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Abstract 

Background: Interferon-beta (IFN-beta) activates the immune response through the type I IFN signaling pathway. 
IFN-beta is important in the response to pathogen infections and is used as a therapy for Multiple Sclerosis. The 
mechanisms of self-regulation and control of this pathway allow precise and environment-dependent response of 
the cells in different conditions. Here we analyzed type I IFN signaling in response to IFN-beta in the macrophage 
cell line RAW 264.7 by RT-PCR, ELISA and xMAP assays. The experimental results were interpreted by means of a 
theoretical model of the pathway. 

Results: Phosphorylation of the STAT1 protein (pSTATI) and mRNA levels of the pSTATI inhibitor SOCS1 displayed 
an attenuated oscillatory behavior after IFN-beta activation. In turn, mRNA levels of the interferon regulatory factor 
IRF1 grew rapidly in the first 50-90 minutes after stimulation until a maximum value, and started to decrease slowly 
around 200-250 min. The analysis of our kinetic model identified a significant role of the negative feedback from 
SOCS1 in driving the observed damped oscillatory dynamics, and of the positive feedback from IRF1 in increasing 
STAT1 basal levels. Our study shows that the system works as a biological damped relaxation oscillator based on a 
phosphorylation-dephosphorylation network centered on STAT1. Moreover, a bifurcation analysis identified 
translocation of pSTATI dimers to the nucleus as a critical step for regulating the dynamics of type I IFN pathway in 
the first steps, which may be important in defining the response to IFN-beta therapy. 

Conclusions: The immunomodulatory effect of IFN-beta signaling in macrophages takes the form of transient 
oscillatory dynamics of the JAK-STAT pathway, whose specific relaxation properties determine the lifetime of the 
cellular response to the cytokine. 

Keywords: Type I interferon pathway, Interferon-beta, Ordinary differential equation, Oscillations, Multiple sclerosis, 
Immunotherapy 



Background 

Type I interferons, such as interferon alfa and beta, are 
cytokines that represent a first-line endogenous defense 
mechanism in response to viruses and bacterial infec- 
tions, are secreted by many cell types (e.g. lymphocytes, 
macrophages and endothelial cells) and they are used as 
a therapy in Multiple Sclerosis (MS). 

The canonical type I interferon (IFN) pathway involves 
different signaling cascades, one of which is the JAK/ 
STAT pathway. This pathway is composed by several 
steps, which include receptor binding (IFNR1 and 2), 
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transformation of the latent transcription factor (a pro- 
tein of the STAT family) into its active form by phos- 
phorylation, nuclear migration of the transcription factor 
(TF), binding of the TF to target promoters, and expres- 
sion of their corresponding genes [1] (Figure 1). Previous 
studies have shown that phosphorylated STAT1 forms 
other TF complexes in response to type II interferons, 
the most important of which is a STAT 1 -STAT 1 
homodimer, known as GAF, that binds to IFN Gamma- 
activated sequence (GAS) elements [2]. 

The target genes of the IFN-beta pathway can be di- 
vided into three categories according to the type of acti- 
vating transcription factor: 1) the ISGF3 complex 
activates genes containing an ISRE binding site in their 
promoter (e.g. ISG15, Mxl, OAS1, IRF7). 2) The GAF 
complex activates genes containing a GAS binding site 
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Figure 1 The canonical type I and type II IFN signalling pathways. The plot represents the canonical IFN pathways and the cross-talk 
between them, including the different pSTAT dimers formed after stimulation. 



in their promoter, such as SOCS1 and IRF1 [3,4]. 3) A 
third class of STAT protein complexes activates other 
canonical pathways that exhibit crosstalk with the JAK/ 
STAT pathway (such as PI3K, NFkB, MAPK) [5]. Re- 
cently it was shown that different immune cell subtypes 
respond differently to IFN-beta induction through acti- 
vation of these different types of genes [6]. 

Different proteins regulate STATl phosphorylation. 
Importantly, a negative feedback loop upon STATl acti- 
vation coexists with a positive feedback mechanism. 
First, the phosphorylation of STATl is inhibited by its 
inhibitor SOCS1 [7]. The SOCS1 protein then inhibits 
STATl phosphorylation at the kinase level. Besides this 
negative loop based on SOCS1, STATl is a subject to 
positive regulation via the TF IRF1, whose transcription 
is induced by activated STATl. IRF1 promotes the ex- 
pression of the STATl gene at the transcriptional level. 
Given the existence of these multiple feedback loops, a 
mathematical modeling of the system would help provide 
an understanding of the response to type I IFN-beta. 

Here we analyzed the type-I IFN-beta signaling pathway 
in macrophages, showing that the response of this path- 
way to IFN-beta stimulation takes the form of transient 
oscillations in STATl phosphorylation. We characterized 
and identified the critical elements governing the transient 



dynamics of IFN activation, and examined the influence of 
this dynamical regime in the response to IFN-beta. 

Dynamical models of IFN induction of the JAK/STAT 
signaling pathway based on nonlinear ordinary differen- 
tial equations, have been previously used to study the ef- 
fect of IFN pre-treatment on the response of the 
immune system to virus infection [8,9] and the robust- 
ness of the pathway to noise and parameter fluctuations 
[10], among other problems. Systems biology approaches 
have also been applied to this pathway in order to exam- 
ine its role in certain pathological mechanisms under- 
lying the behavior of cancer cells [11], and its interaction 
with other key signaling pathways [8,12]. Here we com- 
bine our theoretical model with experimental observa- 
tions. Our results show that a combination of positive 
and negative feedback loops, together with the eventual 
degradation of the IFN signal in the medium, leads to a 
transient oscillatory response in several components of 
the pathway. This behavior is consistent with previous 
numerical results found in pure modeling studies [13], 
and goes beyond previous observations that indicate a 
simpler transient response [10,14,15]. We interpret the 
transient oscillatory response of the pathway in terms of 
the potential effectiveness of IFN-beta treatment in MS 
patients. 
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Results 

IFN-beta induces a transient oscillatory activation of the 
STAT1 pathway 

It is well known that microbial and viral infections in- 
duce endogenous IFN-beta release by macrophages as 
part of the immune cell system response. We could ob- 
serve IFN-beta production accompanied with significant 
increases in levels of phosphorylated STAT1 in the mur- 
ine macrophage-like cell line RAW 264.7 stimulated 
with lipopolysaccharide (LPS) endotoxin and also, with 
viral fragments (poly(I:C)) (data not shown). In this 
study we focused on STAT1 signaling by IFN-beta 
stimulation in macrophages by challenging the RAW cell 
line with increasing concentrations of mouse IFN-beta. 

We observed that phosphorylated STAT1 levels in- 
creased rapidly after IFN-beta induction (Figure 2A). 
The increase was significant as soon as 2 min after 
stimulation and reached a maximum at 10-15 min after 
stimulation, followed by a decrease that correlates with a 



substantial increase in the concentration of the SOCS1 
protein (Figure 2C). A second, smaller peak was visible 
at around 180 min, followed by a subsequent decrease 
back to the baseline level after around 360 min. The 
quick decrease of phosphorylated STAT1 levels is in 
agreement with previous studies pointing to the activa- 
tion of the negative feedback loop mediated by SOCS1 
protein, which suppresses the phosphorylation of JAK1 
and TYK2 proteins and prevents the formation of 
STAT-dimers [16]. The total level of STAT1 protein, on 
the other hand, is maintained practically constant until 
around 200 min after stimulation, after which it starts 
to increase slowly until the end of the experiments 
(Figure 2B). STAT1 mRNA levels grew quickly and 
continuously, starting sharply at around 75 min and 
leveling off after 200 min (Figure 3A). This increase in 
the mRNA level of STAT1, as soon as 1 hour after the 
induction of response by IFN-beta, agrees with the in- 
fluence of the positive feedback loop IRF1 - STAT1 
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Figure 2 Activation of type I IFN pathway by IFN-beta in Raw cell line and the corresponding simulations from the ODE model. A-C) 

pSTATI, STAT1 and SOCS1 protein concentration in cell lysates were measured by Luminex, Flow cytometry or ELISA after stimulation with IFN- 
beta (1000 units/ml). The data was normalized to the maximum level. A) pSTATI protein concentration in the Raw cells stimulated with IFN-beta; 
(B) Total STAT1 protein concentration in the Raw cells stimulated with IFN-beta; C) S0CS1 protein concentration in the Raw cells stimulated with 
IFN-beta. We plot here one out of three independent experiments, which were performed in duplicates. D-F) Model simulation and sensitivity 
analysis with ± 20% change for all parameters (shaded areas) (D: pSTATI; E: STAT1; F: SOCS1). 
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Figure 3 Gene expression levels of STAT1, IRF1 and SOCS1 after IFN-beta stimulation and the corresponding simulations from the ODE 
model. A-C) STAT1, IRF1 and S0CS1 RNA protein concentration in cell lysates were measured by RT-PCR in the Raw cells stimulated with IFN-beta 
(1000 units/ml); A) STAT1 RNA levels in the Raw cells stimulated with IFN-beta; B) IRF1 RNA levels in the Raw cells stimulated with IFN-beta; C) SOCS1 
RNA levels in the Raw cells stimulated with IFN-beta. We plot here one out of three independent experiments, which were performed in duplicates. 
D-F) Model simulation and sensitivity analysis with ± 20% change for all parameters (shaded areas) for STAT1 (D), IRF1 (E), and SOCS1 (F) RNA levels. 



[17] and confirms the importance of this circuit for the 
pathway dynamics. 

To analyze the expression of regulatory genes of the 
type I IFN pathway, we measured the levels of two 
downstream STAT1 genes, SOCS1 (responsible of the 
negative feedback) and IRF1 (mediator of the positive 
feedback). We observed an oscillatory behavior of 
SOCS1 mRNA during the first 360 minutes after stimu- 
lation, with clear peaks at around 90 min and 250 min, 
before returning to baseline levels (Figure 3C). On the 
other hand, IRF1 shows different dynamics, with its con- 
centration raising quickly between 30 and 120 min, then 
reaching a plateau and decreasing more slowly after 
250 min (Figure 3B). We also quantified the expression 
levels of other downstream effector IFN-induced genes, 
such as MX1 and OAS lb, but did not identify any acti- 
vation of their transcription in the RAW 264.7 cell line 
after IFN-beta stimulation (data not shown). These ob- 
servations are in agreement with a differential signal 
transduction mechanism in macrophages when com- 
pared to canonical JAK-STAT pathway in lymphocytes 
[4,18]. Our observations show that in the RAW 264.7 



cell line the main activated genes were the ones con- 
trolled by the STAT1-STAT1 homodimer (IRF1 and 
SOCS1) and containing GAS elements in their promoter 
region. These genes are mainly responsible for the anti- 
microbial activity of the cells [19]. 

Modeling the oscillatory signaling of type I IFN pathway 

In order to identify the critical elements of the signaling 
pathway that are responsible for the transient oscillatory 
behavior observed experimentally in macrophages upon 
IFN-beta stimulation, we built an ODE model based in 
biological knowledge and experimental data (Figure 4). 
To minimize its complexity, we pursued the minimal 
system explaining the experimental observations, instead 
of a full descriptive system [14]. The model was not used 
to reproduce sustained endogenous IFN activation after 
viral infection, although it could be applied to that 
scenario. 

The first ingredient of our model is the binding of 
IFN-beta to the receptor (with the concentration of acti- 
vated receptor being represented by the variable S 
below). The activated receptor induces phosphorylation 
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Figure 4 Mathematical model of type I IFN signaling pathway. 

Graphical representation of the biological processes modeled. 



of the STAT1 protein (represented by A). Phosphory- 
lated STAT1 translocates from the cytoplasm (A pc ) to 
the nucleus (A pn ) (Figure 4). In the nucleus, pSTATl 
complexes activate the transcription of SOCS1 and IRF1 
genes. SOCS1 mRNA (r) is translated into SOCS1 pro- 
tein (R), which inhibits further phosphorylation of 
STAT1. IRF1 mRNA (f) is translated to IRF1 protein (F), 
which activates the transcription of the STAT1 gene 
(with STAT1 mRNA being denoted by a). Each of the 
species has a certain linear degradation rate. We also in- 
troduced receptor internalization through an effective 
degradation (or deactivation) term, consistent with the 
literature [20]. With those ingredients, the model reads: 

dS/dt = b s -\ s S 

dAj dt = b ex pAp n + bdephApc 

b ph SA/k A 



+ b A a- 



dApQ I dt '■ 



l+A/k A + (R/kj) q 
b ph SA/k A 



l+A/kA + iR/h) 
-^statA P c 



-XstatA 

b imp bd e pfoAp C 



dAp n l dt — biynpApc b ex pAp n X-statAp 

A I A, U ( A Pn/ k r)" \ 

dr/dt=br i + (A pn/kr r Xrr 

dR/dt = b R r-\ R R 

{A pn /k f ) n 



df /dt = bj 



f i + (A pn /i< f y 

dF/dt = bpf—\pF 



-V/ 



da/dt = b t 



(F/k F 



u ~^a a + BsTAT 



" 1 + (F/k F ) 

The parameters correspond to transcription rates, in- 
cluding basal transcription of STAT (b a , b p bp B STAT ), 



translation rates (b A , b R , b F ), cooperativity indexes (Hill 
coefficients, n, m, q, u), degradation rates (A/), receptor 
activation (b s ) and deactivation (X s ) rates, phosphoryl- 
ation and dephosphorylation rates (b ph b deph ), and 
nucleo-cytoplasmic transport rates (b imp) b exp ). We ad- 
justed some of the parameters using published data 
sources (Table 1), and the rest were estimated by manual 
fit of the model dynamics to the experimental data 
(Table 2) [21]. Initial conditions are listed in Table 3. 
The model implements the SOCS1 -mediated negative 
feedback loop on pSTATl by means of a competitive in- 
hibition term in the expression determining the phos- 
phorylation rate of STAT1 in the equation for A pc , with 
the parameter kj quantifying the half-maximal inhibition 
threshold and the exponent q defining the sharpness of 
the inhibition. Similarly, the positive feedback through 
IRF1 is described by the transcription activation term of 
STAT1 mRNA in the equation for a, with k F representing 
the half-maximal activation threshold. 

The right panels in Figures 2 and 3 show simulation 
results corresponding to the experimental observations 
presented in the left panels. For comparative purposes, 
both the experimental and model variables were shown 
in relative concentrations dividing by their maximum 
value along the time series. We also performed a sensi- 
tivity analysis by simulating changes of ±20% for every 
model parameter, leading to deviations from the basal 
curve falling within the shaded areas shown in Figures 2 
and 3. 

The model simulations reproduce the main features ob- 
served experimentally, such as the first and very fast peak 
of phosphorylated STAT1 shortly after IFN-beta stimula- 
tion, and the second peak of smaller amplitude at around 
200 min (Figure 2D). The concentration of total STAT1 
protein is almost constant during the first 200 min, after 
which the protein level increases slowly (Figure 2E), fol- 
lowing the increased expression level of STAT1 mRNA, 
which mRNA grows slowly starting at around 50 min after 
stimulation (Figure 3D). In turn, SOCS1 mRNA levels 

Table 1 Parameter values obtained from the literature 



Parameter 


Value 


Reference 


bdeph 


15 min half-life 


[36,37] 


n 


3 


[38] 


U 


1 


[39] 


K 


2.82 hour half-life 


This work 




1.23 hour half-life 


[40] 


A initial 


10 5 molecules/cell 


[41] 


S initial 


700-900 receptors/cell 


[42] 


^SW protein 


24 hour half-life 


[43] 




30 min half-life 


[44] 


f initial 


1 molecule/cell 


[45] 
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Table 2 Parameters for type I IFN ODE model 



Table 3 Initial conditions for type I IFN model simulations 



Name 


Symbol 


Value 


Unity 


Translation rate for STAT1 


b A 


65 


min" 1 


Receptor production rate 


bs 


0 


min" 1 


Basal STAT1 RNA 


BsTATl 


0.0062 


min" 1 


Transcription rate for STAT1 


b a 


0.1 


min" 1 


Transcription rate for S0CS1 


b r 


12.8 


min" 1 


Transcription rate for IRF1 


bf 


2.7 


min" 1 


Translation rate for IRF1 


Of 


I .U ■ I U 


min 


Translation rate for S0CS1 


b R 


1.0- 10 2 


min" 1 


Phosphorylation STAT1 rate 


b ph 


1.3- 10 3 


min" 1 


Dephosphorylation STAT1 rate 


bdeph 


0.036 


min" 1 


Import to the nucleus rate (pSTATI) 


bjmp 


0.013 


min" 1 


Export from the nucleus rate (STAT1) 


bexp 


0.048 


min" 1 


STAT1 phosphorylation activation (Hill's 
constant; half maximal activation) 


k A 


4680 


molecules 


Dissociation constant for the enzyme- 

i « i •_, i r r~\ r~~ r ^\ (\ 1*11/ ■ , i ir 

inhibitor by SOCS1 (Hills constant; half 
maximal activation) 


k. 


82680 


molecules 


SOCS1 transcription activation by nuclear 
pSTATI (Hill's constant; half maximal 
activation) 


1, 


23400 


molecules 


IRF1 transcriptional activation by pSTATI 


k f 


7366 


molecules 


STAT1 transcriptional activation by IRF1 


kp 


1 .3 ■ 1 0 5 


molecules 


Cooperativity of SOCS1 protein over 
STAT1 dimers 


Q 


4 




Cooperativity of STAT1 on SOCS1 gene 
promoter 


n 


3 




Cooperativity of STAT1 on IRF1 gene 

pi Ul 1 lULcl 


m 


2 




Cooperativity of IRF1 on STAT1 gene 
promoter 


u 


l 




Receptor internalization/degradation rate 


As 


0.0229 


min" 1 


SOCS1 RNA degradation rate 


K 


0.0347 


min" 1 


SOCS1 protein degradation rate 


A« 


0.0231 


min-1 


IRF1 RNA degradation rate 


*r 


0.0173 


min" 1 


IRF1 protein degradation rate 


A F 


0.0116 


min" 1 


STAT1 RNA degradation rate 


K 


0.0058 


min" 1 


STAT1 protein degradation rate 


^STAT 


0.0007 


min" 1 



increase from the beginning, showing a first peak at 
around 90 min and a second smaller peak following the 
second peak of phosphorylated STAT1 at around 250 min 
(Figure 3F), in agreement with the experiments. Again 
similarly to the experiments, IRF1 mRNA levels show a 
bell-shaped time course (Figure 3E), with an increase re- 
sembling that of SOCS1 mRNA levels (Figure 3F) and 
remaining high from around 90 min to 250 min, when 
IRF1 mRNA levels decrease to half their maximum value 
at around 350 min. 



Name 


Symbol 


Value 


Non-phosphorylated STAT1 


A 


1.0- 10 5 molecules 


IFN artivatinn rprpntnr 

II IN OLLI VOUUI 1 1 l^yl 


s 


1 0DD mnlpn iIp^ 


Phosphorylated nuclear STAT1 


Apn 


1 molecule 


Phosphorylated cytoplasmic STAT1 


Ape 


10 molecules 


STAT1 mRNA 


a 


1 molecule 


IRF1 mRNA 


f 


1 molecule 


SOCS1 mRNA 


r 


1 molecule 


IRF1 protein 


F 


1 molecules 


SOCS1 protein 


R 


1 molecule 



The model allows us to interpret the second peak ob- 
served experimentally in pSTATI and SOCSl-mRNA 
levels in terms of an underlying damped oscillatory dy- 
namics. We now ask what are the mechanisms leading 
to oscillations, on the one hand, and to damping, on the 
other hand. A well-known gene circuit architecture that 
leads to oscillatory behavior is a combination of positive 
and negative feedback loops [22]. As mentioned above, 
our model contains a negative feedback loop mediated 
by SOCS1. We can examine in the model the effect of 
not having this feedback by eliminating SOCS1 signaling 
from the model. The results show that this negative 
feedback is required for the oscillatory behavior to arise: 
its absence leads to a transient plateau of high pSTATI 
levels during the first 4-5 h of IFN treatment (Figure 5B), 
which contrasts with the relaxation oscillator behavior 
obtained for our basal parameter values (Figure 5A), 
which is a closer match to the experimental observations 
(Figures 2-3). The model also contains a positive feed- 
back loop mediated by IRF1. This loop, however, does 
not appear to be crucial for the oscillatory behavior of 
pSTATI (Figure 5C), and is only necessary to reproduce 
the experimentally observed increase of STAT1 expres- 
sion (Figure 3A D). 

The combination of negative and positive feedbacks 
discussed in the preceding paragraph would naturally 
lead to sustained oscillatory behavior. The experimental 
observations shown in Figures 2 and 3, however, reveal a 
strong damping of the oscillations that lead to their sud- 
den disappearance. This behavior is not consistent with 
the standard damping undergone by nonlinear oscilla- 
tions when they become unstable via a Hopf bifurcation, 
in which case the damping is either slow close to the bi- 
furcation, or the damped oscillations are too weak to 
begin with far away from the bifurcation. A key distinct- 
ive characteristic of our model is the fact that the exter- 
nal input to which the system is subject (mediated by 
the activated receptors represented by S in the model 
above) decays monotonously due to receptor inactiva- 
tion by internalization or degradation [20,23]. Assuming 
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Figure 5 Type I IFN pathway model simulations during the first 8 hours after IFN-beta stimulation. (A) Model simulations showing 
oscillations of total (nuclear plus cytoplasmic) pSTATI protein, SOCS1 mRNA expression (dashed line) and STAT1 mRNA expression (dotted line). 
(B) Corresponding simulations where the SOCS1 -mediated negative feedback is disrupted by assuming an infinite value of the repression 
threshold k,. (C) Corresponding simulations not including the IRF1 -mediated positive feedback, by assuming a zero value of the STAT1 activation 
threshold k F (which leads to saturation of the corresponding Hill function, so that the dependence of 0 expression on F is removed). 



a linear decay, the external input decreases exponentially 
(Figure 6A), taking the system quickly out of the oscil- 
latory regime and thus leading to a strong damping of 
the oscillations, as seen experimentally. Exploring sys- 
tematically the dependence of the dynamics on the re- 
ceptor level, we observed that as the initial levels of 
activated receptor decrease (Figure 6A) the second peak of 
both pSTATI (Figure 6B) and SOCS1 mRNA (Figure 6C) 
levels diminish, with the SOCS1 expression peak dis- 
appearing earlier than the pSTATI concentration peak. 

Bifurcation analysis of the ST ATI pathway model 
identifies translocation to the nucleus as a critical step 

The temporal evolution of the phosphorylation of 
STAT1 can be crucial for understanding the response to 
IFN-beta therapy, and may provide an explanation of the 
lack of response to this therapy in some cases [24]. In 
particular, transient oscillatory dynamics could provide a 
way for the STAT1 pathway to increase the duration 



of its response to IFN-beta in a physiological manner 
(i.e. without a period of sustained constant activation 
as in Figure 5B). In order to establish the conditions 
under which this transient dynamics exists, we ana- 
lyzed the behavior of the system for combinations of 
two-parameter pairs, distinguishing between the param- 
eter values for which pSTAT dynamics is overdamped 
(and thus non-oscillatory) and those for which the oscilla- 
tions are underdamped (which corresponds to the experi- 
mental situation reported above). We focused on the 
phosphorylation and dephosphorylation rates (b ph and 
bdeph) and export and import rates (b exp and b imp ). These 
parameters represent crucial steps to regulate the nuclear 
availability of transcription factors such as pSTAT dimers, 
and thus also the expression of downstream genes. 

In order to identify mathematically the two different 
dynamical regimes mentioned above, we examined the 
stability properties of the steady state of the ODE system 
for a constant activated receptor level (S) equal to its 



Pertsovskaya et al. BMC Systems Biology 201 3, 7:59 
http://www.biomedcentral.eom/1752-0509/7/59 



Page 8 of 12 




Receptor activation (S) 



200 300 

time [min] 



400 



500 



B 



o 

1 

|i 

o 



I ? 



pSTATI protein 




200 300 

time [min] 



c 
c 



£ § 

* c 
c 



00 



1 ? 



S0CS1 mRNA 



~ 08 



06 



04 



0) 0.2 



S-'OK 



°0 



100 



200 

time [min] 



:::::;:Y***^ 

300 400 



500 



Figure 6 Influence of activated receptor level on the transient oscillatory dynamics. (A) Time evolution of the activation level of type I IFN 
receptors in the model (S) by stimulation with IFN-beta (added at t = 0) for different initial levels (in red, S = 500 molecules at t = 0). (B, C) pSTATI 
(B) and SOCS1 mRNA (C) dynamic responses for varying levels of initial receptor activation as in panel A (lower lines in B and C correspond to 
lower lines in A). 



initial value. In that context, the underdamped/oscillatory 
regime is characterized by a steady state that takes the 
form of an unstable focus (i.e. the stability eigenvalue with 
maximum real part has negative imaginary part), whereas 
the steady state in the overdamped/ non-oscillatory regime 
is a node (the stability eigenvalue with maximum real part 
has no imaginary component). In that way, by calculating 
the imaginary part of the stability eigenvalue of the steady 
state with maximum real part, we can identify the param- 
eter regions in which the pathway exhibits a transient os- 
cillatory response to IFN-beta. The result, for the 
parameter space formed by the phosphorylation and de- 
phosphorylation rates b p h and bd ep h> is shown in Figure 7 A. 
This figure shows, on the one hand, the prevalence of os- 
cillations for a wide range of these parameters, and on the 
other hand it tells us the conditions for which transient 
oscillations exist. For instance, increasing sufficiently the 
dephosphorylation rate can transform an oscillatory 



regime into a non-oscillatory one, and reversely, by mak- 
ing the phosphorylation rate large enough the system can 
be made to exhibit transient oscillations. Figure 7B shows 
examples of these two dynamical regimes for two specific 
parameter sets within this phase diagram. 

We also tested the influence of the import rate of 
pSTATI molecules into the nucleus, and of the export 
rate of STAT1 from the nucleus into the cytosol. By tun- 
ing both parameters (b exp and b imp ) simultaneously, we 
observed again that the transient oscillatory regime is 
prevalent in this system (Figure 7C). We found that the 
oscillatory regime is associated with high nuclear import 
rates in combination with high export rates. For high ex- 
port rates but low import rates, the pathway exhibits an 
overdamped (non-oscillatory) response, showing a 
sustained plateau in the transient level of pSTATI (dis- 
continuous blue line in Figure 7D). Conversely, for high 
import rates and low export rates the response is also 
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time (min) time (min) 

Figure 7 Stability analysis of the steady state solution in two 2D parameter spaces. (A, C) The color scale represents the absolute value of 
the imaginary part of the stability eigenvalue with maximum real part, corresponding to the steady state of the system after IFN stimulation, for 
varying phospho/ dephosphorylation rates {b ph and b dephl panel A), and nuclear export/import rates {b exp and b impi panel B). Two distinct 
dynamic regimes can be identified with this analysis: the damped oscillatory regime (shifted to red) and the overdamped/ non-oscillatory regime 
(blue). (B, D) Examples of pSTATI time evolution in both regimes (damped and over-damped in red/blue lines, respectively) corresponding to 
parameter position of circle markers for the two diagrams shown in panels A and B, respectively. 

V . J 



overdamped, but with a faster decay (continuous blue 
line in Figure 7D). 

Discussion 

The aim of this work was to characterize the dynamics 
of the key components of the type I IFN-beta signaling 
pathway in macrophage RAW 264.7 cells. This system 
robustly translates extracellular chemical signals through 
cell membrane receptors, leading to phosphorylation of 
the STAT transcription factors, which induce gene ex- 
pression of multiple targets. JAK/STAT signaling directly 
regulates the immune system response under viral or 
bacterial infection, and is also important in auto- 
immune diseases and cancer treatments. The IFN sig- 
naling network affects different complex pathways, 
involving processes such as differentiation, proliferation, 
survival and cell death. Importantly, it is a canonical 



pathway involved in first-line treatments of multiple scler- 
osis as a main target of the IFN system [25] but, also, af- 
fects different complex pathways, involving processes 
such as differentiation, proliferation or survival and 
cell death [26,27]. 

In this paper, we used a combination of experimental 
approaches in order to obtain a quantitative picture of the 
response of the JAK/STAT signaling pathway to IFN-beta 
stimulation, and to identify the most relevant aspects of 
its dynamics to be modeled with kinetic equations. Experi- 
ments uncovered several important features of JAK/STAT 
signaling dynamics during the first eight hours after treat- 
ment with IFN-beta. For example, our results showed the 
transient oscillatory nature of STAT1 activation (pSTATI), 
with a fast increase in cytosol concentration early after 
stimulation (within the first hour), followed by a second- 
ary concentration peak at around 200 min. A key STAT1 
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transcription target such as SOCS1 also showed two peaks 
of expression (correlated in time to the pSTAT peaks) at 
around 90 min and 250 min after stimulation, whereas an- 
other important target, namely IRF1, exhibited a more 
bell-shaped plateau signal, respectively (Figures 2 and 3). 
Our model simulations also exhibit a transient oscillatory 
behavior in pSTATl concentration, and reveal that the os- 
cillations require the presence of a negative feedback loop 
on STAT1, mediated by its phosphorylation inhibitor 
SOCSL Previous mathematical models of the type I and 
type II IFN pathways have suggested the possibility that 
STAT1 pathway has an oscillatory behavior [9,13] and in- 
dicated the importance of the SOCS1 negative feedback 
[10,14,28,29]. Another factor that has been proposed to be 
important in defining the response to IFN is the basal 
level of receptors of the JAK/STAT pathway [30]. In our 
model this aspect was also taken into account, showing 
clear effects on the dynamics of the pathway response 
(Figure 6). Going beyond previous models, our theoretical 
results show that the physiological regime of the pathways 
response to IFN-beta takes the form of damped oscilla- 
tions that can be identified by means of a stability analysis 
of the models steady state solution. This analysis shows 
that processes such as the phosphorylation and dephos- 
phorylation of pSTATl, and the transport of STAT1 be- 
tween the nuclear and cytosol compartments, can make 
the pathway switch between underdamped and 
overdamped oscillatory regimes [31]. 

Implications of the type I IFN signaling dynamics in 
autoimmune diseases 

IFN-beta is the most common treatment for MS [32], 
exerting a pleiotropic immunomodulatory activity not 
well understood [25]. IFN-beta treatment decreases acti- 
vation, proliferation, cytokine release, and migratory 
properties of activated T cells, diminishing their ability 
to enter and damage the brain tissue. In spite of these 
properties, up to 40% of patients do not respond to IFN- 
beta therapy, which represents a significant health prob- 
lem [24]. Previous genomic studies have identified certain 
genes belonging to the IFN pathway that are associated 
with a lack of response to IFN-beta, suggesting that the 
genetic background of certain individuals may modulate 
this pathway, and consequently the response to therapy, 
by specific transcriptional profiles [33]. For example, it 
was recently shown that the response to IFN-beta differs 
between immune cells, and an analysis of non-responders 
to IFN-beta therapy indicates an impairment of the type I 
IFN pathway in the monocytes of those patients [6,34]. 

Our study indicates the importance of identifying the 
temporal dynamics of the concentration of certain key 
components of the JAK-STAT pathway, such as the 
phosphorylated form of the STAT1 protein, and of the 
expression of interferon-stimulated transcription genes 



like SOCS1 and IRF1, within the first 8 hours of IFN- 
beta administration. Cataloguing these dynamics could 
provide us with early molecular biomarkers that allow us 
to distinguish the lack of response to IFN-beta therapy 
of certain MS patients. 

Methods 

Materials and reagents 

Cells were obtained from ATCC library, mouse recom- 
binant IFN-beta was purchased from Cell sciences, lipo- 
polysaccharide from Escherichia coli and poly (I: C) salt 
was purchased from Sigma- Albrich, lipofectamine 2000, 
Hiperfect transfection agent were purchased from 
QIAGEN, Taqman PCR master mix, VIC-dye GAPDH 
endogenous control, IRF1, SOCS1, STAT1, STAT2, 
MX1, OAS la pre-designed FAM-dye assays were pur- 
chased from Applied Biosystems, total STAT1 and 
STATl(pTyr701) antibodies and beads, cell detection kit 
for xMAp assays were purchased from Merck Millipore 
(Billerica). Alexa Flour 647 STAT1 (pTyr701) and PE 
STAT1 N- terminal anti-Mouse antibodies and all buffers 
for cytometry were purchased from BD biosciences. 
APC-labelled IFNAR1 antibody was purchased from 
Biolegend. 

Cell culture and stimulation 

Mouse leukemic monocyte macrophage cell line RAW 
264.7 cell line was purchased from ATCC and maintained 
in DMEM medium complemented with 10% fetal bovine 
serum and 1% antibiotics at 37°C and 5% C02. The cells 
were passed every 2-3 days and maintained in 20-80% 
surface coverage. One day before the stimulation the cells 
were seeded in 12 well plates in concentration 1 x 10 6 
cells/ well. The cells were stimulated with 1000 units of re- 
combinant mouse IFN-beta, 15 \i% of LPS for different 
times or 15 (ig of poly (I:C) solution. At the end of stimu- 
lation supernatants or cells were collected for further ana- 
lysis. The same amount of PBS was added at the 
corresponding time-points to the control samples. 

RT-PCR 

Cell lysates were prepared with QiaShredder columns 
and total RNA was isolated using standard Qiagen 
Rnaesy Mini kit protocol. Equal amount of total RNA 
was added to each reverse transcription reaction tube 
(High-Capacity cDNA Reverse Transcription Kit from 
Applied Biosystems and cDNA was used for a second 
step of RT-PCR. Results were analyzed using relative 
2CTT method normalized to a GAPDH endogenous 
control (VIC-dye primer-limited control from Applied 
Biosystems) as described before [35]. All the qRT-PCR 
experiments were performed in triplicates and repeated 
three times independently. 
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Western blot and quantification 

Western blot (WB) was performed using polyclonal 
rabbit anti-mouse pSTATl and STAT1 N- terminal anti- 
bodies (Abeam) using standard WB protocol Western 
blot results were quantified using ImageJ software 
(http://rsb.info.nih.gov/ij/index.html) using the method 
of Luke Miller (http://www.lukemiller.org/journal/2007/ 
08/quantifying-western-blots-without.html) 

ELISA and xMAP multiplexing assays 

IFN-beta in culture supernatants and SOCS1 protein 
concentration in cell lysates were measured by standard 
sandwich ELISA with anti-mouse SOCS1 antibodies 
(Abeam). IRF1 protein concentration in cell lysates was 
measured by in-cell ELISA using the kit (Thermo Scien- 
tific) STAT1 total protein and phosphorylated state 
(Tyr701) concentrations (nuclear and cytoplasmic to- 
gether) were measured using xMAP assays and read in 
Luminex 201 platform using standard vacuum separ- 
ation protocol (Millipore). xMAP experiments were re- 
peated twice. 

Flow cytometry 

Cells for flow cytometry were stimulated with 1,000 un/ml 
of IFN-beta as stated before and fixed immediately after 
stimulation. IFNAR1 receptor on the surface of the RAW 
264.7 cells was marked using APC-labelled anti-IFNAR anti- 
body. The mean fluorescent intensity was calculated using 
Flowjo software. For STAT1 staining cells were fixed imme- 
diately after stimulation in order to preserve phosphoryl- 
ation and then permeabilized using Perm III buffer (BD 
biosciences). Samples were stained simultaneously with anti- 
STAT1 (pTyr701) and anti-STATl total (N-terminus) anti- 
bodies. The mean fluorescent intensity, the percent of 
staining-positive cells, the medians and the standard devi- 
ation were calculated using Flowjo software and the raw 
single-cell data were extracted to plot the histograms and 
further analysis. 

Mathematical modeling 

The model and simulations were run in MATLAB using the 
ODE15s solver (Matlab codes are provided in the Additional 
files 1, 2 and 3). The stability analysis of the dynamical sys- 
tem was performed with custom-made Matlab codes. 

Availability of supporting data 

The model is available as a matlab script in the supporting 
materials. The raw experimental data are available from 
the authors upon request. 

Additional files 



Additional file 1: Type I IFN pathway ODE model equations. 



Additional file 2: Type I IFN pathway ODE model as Matlab file. 

Additional file 3: Integration and solutions of the Type I IFN 
pathway ODE model as Matlab file. 
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